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CHAPTER  I 

if!  , 

Introduction 

-i/u  J 

In  this  paper  we  shall  develop  two  algorithms  for  curve  fitting 
which  will  be  adaptive  in  nature  using  results  from  uniform  approxi- 
mation theory.  The  organization  of  this  paper  is  as  follows:  in  this 

section  we  shall  review  some  basic  results  from  uniform  approximation 
theory,  uniform  approximation  theory  with  side  conditions,  and  the 
Remes  algorithm;  in  the  second  section  we  shall  give  our  curve  fitting 
algorithms;  and  in  the  last  section  we  shall  report  on  our  numerical 
testing  of  our  algorithms. 

Let  X be  a compact  topological  space.  Let  C(X)  denote  the  family 
of  all  real-valued  continuous  functions  defined  on  X.  As  is  well 
known,  C(X)  is  a complete  normed  linear  space  under  the  supremum 
(Chebyshev,  uniform,  max,  infinity)  norm 

||f||  = max  |f  (x)  | , f L C(X). 
xf  X 

Let  S be  a subset  of  C(X)  and  f € C(X)  be  fixed.  Consider 

inf  ||f  - p ||  = d. 
p«S 

If  there  is  a p*  € S such  that  ||p*  - f ||  = d,  then  p*  is  called  a best 
approximation  to  f from  S.  For  a given  subset  S,  the  following  ques- 
tions may  be  posed. 

1.  Does  a best  approximation  exist  for  a particular  f € C(X)? 

For  each  f € C(X)? 
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2.  If  a best  approximation  does  exist,  is  it  unique? 

3.  If  a best  approximation  exists,  is  there  an  efficient  way  of 
computing  it? 

4.  Does  a best  approximation  depend  continuously  on  the  function 
being  approximated?  (That  is,  will  small  errors  inevitably 
introduced  when  calculating  best  approximations  prevent  an 
algorithm  from  succeeding  in  converging  to  something  close  to 
a best  approximation?) 

In  order  to  answer  these  questions,  we  restrict  our  attention  to 
the  following  X and  S.  Let  X be  a compact  subset  of  an  interval  on  the 
real  line,  and  S be  a (finite  dimensional)  Haar  subspace  of  C(X).  That 
is,  if  S has  dimension  n,  then  the  only  function  in  S having  n or  more 
zeros  is  the  zero  function. 

For  this  particular  X and  S,  all  of  the  above  questions  have  been 
answered  in  the  affirmative.  For  theorems  and  proofs  regarding  ques- 
tions 1,  2 and  4,  see  for  example  [33.  We  now  address  ourselves 
specifically  to  question  3.  The  following  theorem  characterizing  best 
uniform  approximations  will  prove  useful. 


CHARACTERIZATION  THEOREM  [3,  p.  753*  If  S has  dimension  n,  then 

p*  C S is  a best  approximation  to  f € C(X)  ^ S if  and  only  if  there  are 
n + 1 points,  x^  ....  x^+1  in  X C [a,  b3  with  x.  < xi+1»  (i=l,...,  n) 
such  that  f(x£)  - p(Xj)  = (-l)1*,  and  |x|  = ||f  - p*||  . 

S I Given  any  f G C(X)  ^ S,  then  any  set  of  n + 1 distinct  points  from 


X determines  uniquely  a p € S and  a X € R so  that  p(x^)  + (-D^X 
= f(x.)  since  S is  a Haar  subspace  (see  [33).  The  characterization 


theorem  and  the  above  observation  allow  us  to  view  the  problem  of 

finding  the  best  approximation  to  f £ C(X)  as  finding  a set  (called  an 

extreme  or  alternating  set)  of  n + 1 points,  (x.  , ...»  x , }C  X where 

t n+1 

xi  < x^+1>  i = 1»  •••»  n so  that  the  p and  the  A determined  by  the 
extreme  set  have  the  property  that  |x|  = ||f  - p||  . In  this  case,  p 

must  be  the  best  approximation  to  f from  S. 

We  now  describe  Remes  second  (or  multiple  exchange)  algorithm 
which  calculates  the  best  approximation  to  f € C(X)  'v  S by  finding  an 
extreme  set  for  f.  We  assume  that  f is  fixed  and  that  an  arbitrary 
reference  set  R = (x^,  ...,  x^+^}  of  n + 1 distinct  points  from  X has 
been  chosen  so  that  x^  < x^+1»  i = 1.  •••»  n,  with  the  added  assumption 
that  the  X determined  by  this  reference  set  is  not  zero.  Since  f ^ S, 
and  S has  only  dimension  n,  such  a set  may  always  be  readily  found. 

We  begin  each  iteration  by  computing  the  p and  X associated  with 
the  reference  set  R.  Let  o(x)  = f(x)  - p(x).  Since  a(x.)  = (-1)1X 
with  A * 0,  a(x)  changes  sign  from  one  extreme  point  to  the  next.  If 
[x^,  x.+^]  C X,  then  since  f and  p are  both  continuous  a has  a zero 
in  Cx^ , Xi+1L  call  it  z^.  If  [x^  xi+1]  £ X,  let  z.=min{x£X:  x > xi 

and  o(x.)*o(x)  <_  0).  Set  zQ  = min{x  £ X),  zn+1  = max{<  € X}.  (It  is 
possible  that  zn  = zn+^,  although  this  will  present  no  difficulty  as 
will  become  clear  later). 

Now  we  construct  a new  reference  set  Y = {y^  ....  yn+1)  C X of 

distinct  points  so  that  the  function  a(x)  still  alternates  in  sign, 

that  is  o(yj)*o(y^+1 ) < 0 for  i = 1,  ...»  n,  |o(y^)|  >_  |o(x.)|  for 

i = 1,  . n ♦ 1,  and  for  some  k,  k = 1,  . ...  n + 1 |o(yk)|=  ||  f -p  ||  . 

We  first  find  a trial  reference  set  ^ = (y^,  ....  Yn+1^  by  choosing  y\ 
to  be  a point  in  [z^,  z^]  f)  X where  o(x)*(sign(o(xi)))  attains  its 


4 


maximum  on  [z^  z^]  fl  X.  If  there  is  a k such  that  | a(y^ ) | = ||f-p||  « 

then  we  set  Y = Otherwise,  let  y*  be  a point  in  X so  that  | o ( y * ) | 

= 11^  ~ P ||  • We  form  Y be  adding  y*  to  If  and  removing  one  of  the  y^'s  so 
that  Y is  still  ordered  and  a(x)  alternates  in  sign  on  Y.  This  may  be 
accomplished  as  follows.  If  y*  £ (y^,  y )0  X,  choose  j to  be  either 
k or  k+1  so  that  o(yj)  is  of  the  same  sign  as  a(yA).  Then  replace  y^ 

f\j 

by  y*  to  form  Y.  If  y*  < y^,  eliminate  yn+^  (by  the  way  y^  was  chosen, 

it  must  be  the  case  that  signCcKy^))  = -sign(o(y  + 1)).  Similarly,  if 
. *\»  . . 'I 

y*  > y , eliminate  y^. 

If  the  new  reference  set  Y is  the  same  as  the  reference  set  from 
the  last  iteration,  R,  then  it  is  necessarily  true  that  | a | = ||f  - p||  , 
and  from  the  Characterization  theorem,  p must  be  the  best  approximation 
to  f from  S,  and  we  terminate  the  algorithm.  If  Y t R,  then  we 
iterate  again,  substituting  our  new  reference  set  Y for  the  old  refer- 
ence set  R.  The  exchange  procedure  outlined  above  insures  that  as  long 
as  the  A determined  by  the  initial  reference  set  is  non-zero  then  at 
each  successive  iteration,  the  magnitude  of  the  present  A is  strictly 
greater  than  the  magnitude  of  the  A in  the  preceding  iteration  (see  [3]). 

It  can  be  shown  [6,  p.  108]  that  the  above  algorithm  does  converge 
(uniformly)  to  the  best  approximation  to  f from  S.  In  fact,  under 
stronger  hypotheses,  the  convergence  is  quadratic  [6,  p.  111]. 

We  may  extend  these  results  by  forcing  the  approximation  functions 
to  satisfy  certain  side  conditions.  We  now  consider  two  such  exten- 
| sions.  First,  we  consider  imposing  interpolatory  constraints  on  the 

approximating  functions  and  their  derivatives,  then  we  shall  consider 
restricting  the  range  of  the  approximating  functions.  In  both  cases, 
if  we  make  strong  enough  assumptions  on  the  structure  of  the  subspace  S, 

\ 
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the  questions  of  existence,  uniqueness  and  continuous  dependence  on 
the  data  can  be  affirmatively  answered  (see  [5]  for  interpolatory 
constraints;  see  [8]  for  restricted  range).  Moreover,  in  both  cases, 
the  Remes  algorithm  can  be  extended  to  effectively  compute  best 
approximations  from  these  new  approximation  classes. 

In  order  to  consider  the  interpolatory  constraint  problem  we  must 
first  define  an  extended  Haar  system.  We  follow  the  notation  in  [5]. 
Let  (u^(t)}?_1  be  a family  of  functions  in  C[a,  b]  (with  each  suffi- 
ciently differentiable  so  that  what  follows  is  well  defined)  and 
[a,  b]  be  such  that  a ^ tg  1.  • • • .1  t ^ b.  Define 


u*(i* • • • »n  \ - 

S. V 


'Y, 

ui ( 1 1 ’ W 


V*l>  W 


Ul(tn) 


2f_(t  ) 
n n 


where  for  fixed  j 


Vv 


...  = t. 
3 


1 < i < n.  The  functions  (u^(t)}?_^  will  be  called  an  extended 
Chebyshev  system  of  order  v on  [a,  b]  provided  u^  £ Cv  ^[a,  b], 
i = 1 , . . . , n , and 


u*(J ) > 0 

T1 n 

for  all  choices  <_  tj  <_  ...  <_  tn>  t . £ [a,  b],  where  equality  occurs 
in  groups  of  at  most  v consecutive  tj  values.  An  n-dimensional  subspace 


I 
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M of  C[a,  b]  is  said  to  be  an  extended  Haar  subspace  of  order  v, 

v n,  provided  that  there  exists  a basis  for  M which  is  an  extended 

Chebyshev  system  of  order  v on  [a,  b]. 

We  may  state  the  interpolatory  constraint  problem  as  follows. 

Let  M be  an  n-dimensional  extended  Haar  subspace  of  order  v (on 

[a,  b]).  Let  c X C [a,  b]  satisfy  x1  < x?  < ...  < xk>  and 

let  (m. }.  be  a set  of  positive  integers  such  that  max  (m.)  - 1 

l<i<k  1 
k 

£ v,  and  £ m.  = m < n.  We  assume  X contains  at  least  n - m + k + 1 
i-1  m. _i 

points.  Let  be  a set  m real  numbers  and  define 

/ • \ 

S = {p  € M : p ^ (x. ) = a. . , 1 <i  <k  and  0 < i < m.  - 1}. 

l — — — — i 

The  interpolating  constraint  problen  then  is  to  find  best  approximations 
to  functions  in  the  class 


2f(x)  = {f  e C(x)  : f(x;.)  = a.Q,  1 <_  i <_  k} 


from  the  class  S. 

The  Remes  algorithm  may  be  extended  in  the  following  manner  to 
compute  interpolatory  approximations.  Let  f € fr(x),  {x. 

(m.  }*  and  {a..}^’"1*  * be  as  above.  Let  R = (t,,  ...»  t _}C  X 
i i=l  1 j 1=1 , j =u  l n+1 

be  such  that  t^  ^ ^ t , and  each  x^  is  repeated  in  R m. 

k 

times.  Moreover,  if  t j € R 'v  U x^,  then  tj_^  < tj  < tj+l‘  ^et 

^ be  such  that  R 'v  U It  can  be  shown  (see 

[5])  that  there  exist  a unique  p € S and  a real  number  X satisfying 


< 
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and 


f(t.  ) - p(t.  ) = (-1)  l\,  1=1,  2,  n-m+1 , 

J*. 

p(^(x.)  = a . ^ , j = 0,  1,  ....  m..  - 1,  i=l,...,k. 


We  assume  that  R is  initially  chosen  so  that  | A j >0. 

Next  we  perform  an  exchange  on  the  reference  set  R to  form  a new 
reference  set  Y C X with  the  properties  that  each  x^  is  repeated  in  Y 

m^  times,  that  there  is  a y*  G Y such  that  | f (y* ) - p(y*)|  = || f — p 1 1 , 

s-  k k 

and  if  (t,  = Y s,  U x.  then 

A J 


sign{f (t . ) - p(t.  )}  = (-l)q  sign{f(t.  ) - p(t.  )},  l - 1,...,  n-m+1, 

. H 11  *1 

where  q = i,  - i^.  After  the  exchange,  if  R ^ Y we  iterate  again  using 
the  reference  set  Y in  place  of  R.  If  R = Y then  the  present  p is  the 
best  approximation  to  f from  S,  and  the  algorithm  is  terminated.  For  a 
description  of  the  exchange  and  a proof  that  this  algorithm  converges 
uniformly,  see  [5]. 

We  now  consider  uniform  restricted  range  approximations.  For 
a given  f € C(X)  we  assume  that  we  have  two  extended  real  valued 
functions  u(x)  and  £(x)  defined  on  X which  satisfy  the  following 


requirements : 

(i) 

l(x) 

may  assume  the  value  but  never  +<» 

(ii) 

u(  X ) 

may  assume  the  value  +®  but  never 

( iii ) 

A = 

l and  B = u 1(®)  are  open  sets  : 

(iv) 

fc(x) 

is  continuous  on  X 'v  A. 

(v) 

u(x ) 

is  continuous  on  X 'v  B. 

(vi) 

*(x) 

£ f(x)  <_  u(x)  for  all  x € X. 

(vii) 

A(x) 

< u(x)  for  all  x £ X. 
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Let  S be  an  n-dimensional  Haar  subspace  of  C(X).  Let 
K = {p  € S : i(x)  p(x)  u(x)  for  every  x €.  X}.  We  say  that  p*  is  a 
best  restricted  range  approximation  to  f if  p*  £ K and  || f-p* ||  =inf  ||f-p||. 


ptK 


Taylor  [8]  has  proved  that  under  the  above  hypotheses  on  f,  £,  u and  S, 
and  if  K i <p,  then  there  exists  a unique  best  restricted  range  approxi- 
mation to  f.  As  long  as  there  is  more  than  a single  element  in  K,  then 
we  have  a characterization  of  the  best  restricted  range  approximation  to 
f as  given  in  the  following  theorem.  Before  starting  Taylor's  charac- 
terization theorem,  we  define  for  each  p £ K the  following  subsets  of  X. 

X+1  = {x  € X : f(x)  - p(x)  = ||  f - p ||  > 


X_1  = {x  £ X : f(x)  - p(x)  = 

X+2  = {x  £ X : p(x)  = JL(x)> 

X_2  = {x  fe  X : p(x)  = u(x)} 

Xp  ■ X.1U  X-1U  X«U  X-2* 


P||) 


Restricted  Range  Characterization  Theorem  [8].  p*€  K C S is  a best 

restricted  range  approximation  to  f if  and  only  if  there  exist  n + 1 
distinct  points,  x^,  ...,  xfl  with  xQ  < x1  < ...  < xn  in  X^ft  satisfying 
a(x^)  = (-l)1+1a(x1),  where  o(x)  = -1  if  x £ X^U  X_2  and  o(x)  = +1  if 
x€  XtiU  Xf2. 

This  characterization  is  the  basis  for  the  following  multiple 
exchange  algorithm  due  to  Gimlin,  Cavin  and  Budge  [4]  which  is  the 
analog  of  the  multiple  exchange  Remes  algorithm  previously  described.  We 
assume  f € C(X)  is  fixed  and  that  u,  l,  and  K have  been  chosen  as 
described  above.  We  also  assume  that  we  have  an  arbitrary  reference  set 
R = {x^,  ....  *n+^}  C X with  x^  < x2  < ...  < on  which  no  function 

in  K interpolates  f.  Set  oi  = (-l)i+1,  ^ = ffx^,  i * 1,  .... 


n + 1. 
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We  begin  the  first  iteration  by  finding  p£  S and  a real  number 

such  that  p(-x^)  + = f^.  We  next  execute  an  exchange  algorithm  which 

returns  a new  reference  set  Y C X,  and  (possibly)  new  values  for  o.  and 
% 

f^.  If  our  current  reference  set  R ^ Y,  we  iterate  again  substituting  Y 
for  R.  If  R = Y,  then  our  present  p is  in  K and  satisfies  the  character- 
ization theorem  so  that  p is  the  best  restricted  range  approximation  to 
f. 

We  now  describe  the  exchange  algorithm.  We  define  the  function 
sign*  as  follows. 


{sign(f(x)-p(x)  if  i(x)  < f(x)  < u(x)  or  f(x)  i p(x)) 
+1  if  f(x)  = £(x)  = p(x) 

-1  if  f(x)  = u(x)  = p(x) 

It  can  be  shown  (see  [4])  that  at  each  iteration  sign*(f(x^)  - p(x^)) 

= -sign*(f(xi+1)  - P(xi+1))>  i = 1»  ....  n.  Set  zi  = inf{x  £ [x^x.^lOX: 
sign*(f (x)-p(x) ; -sign*(f (x- )-p(x^ ) ) £0 },  i = 1,2,. ..,n.  Set  zQ=min{x  € X), 
zr+^  = max(x  €.  X).  Define  E,.,  , hk  as  follows. 

E.  = max  (s.*(f(x)  - p(x))  - | A | ) 

X£tZi_i,Zi]nX 

= max  (Si*(u(x)  - p(x))} 

xe[z^_i  ,Zi  ]nX 

A 

m.  = max  {S.*(i(x)  - p(x)0 

x«[z.  .z^flX 

where  S^=sign*(  f (x . )-p(  x^ ) ) . Let  Y|=max(E^.  ,M.  ,m^  } , in  case  of  equality, 
choose  Y-  to  be  the  first  largest  member  of  the  triple.  We  define  a new  trial 

f\.  /\#  A 

reference  set  i = (y y^}  by  taking  y^  to  be  a point  fej  ^z.JOX  at  which  y . 


occurs. 
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Define  y = max{E,  M,  m}  where 

E = max^ I x)  - p(x)|  - |x|} 
xeX 

M = max{p(x)  - u(x)} , 
xfeX 

m = max{i(x)  - p(x)}. 
xeX 

In  case  of  equality,  choose  y to  be  the  first  largest  member  of  this 
triple.  If  for  some  k,  1 _<  k n,  Yk  = y,  then  set  Y = Y.  Otherwise,  let 
y*  be  a point  at  which  y occurs,  insert  y*  into  Y and  delete  a y^C  Y 
exactly  a was  done  in  the  Remes  algorithm,  preserving  sign*  alternation. 
We  now  determine  our  new  sets  and  If  y - £ Y is  a point 

at  which  E or  Ej>  occurs,  set  ck  = signffCy.^)  - p(yi>),  and  set  f.=f(y^). 

If  y^  is  a point  at  which  M or  Ph  occurs,  set  ck  = 0,  and  set  f^=u(yi>. 

If  y^  is  a point  at  which  m or  m^  occurs,  set  ck  = 0,  and  set  f^  = Jl(yi). 

This  completes  the  exchange  algorithm. 

See  Gimlin  [4]  for  a proof  that  this  algorithm  does  converge  to  the 
best  restricted  range  approximation  to  f. 

As  can  be  done  with  best  uniform  approximations,  best  restricted 
range  approximations  can  be  computed  subject  to  interpolatory  constraints. 
The  precise  formulation  of  the  problem  is  completely  analogous  to  that 
of  best  uniform  approximations  with  interpolatory  constraints.  A Remes- 
like  algorithm  for  computing  restricted  range  approximations  subject  to 
interpolatory  constraints  is  available  and  is  analogous  to  the  Remes- 
like  algorithm  given  above  for  computing  best  uniform  approximations  with 
interpolatory  constraints.  For  a complete  discussion  of  this  problem  see 
[2]. 


CHAPTER  II 


Two  Adaptive  Piecewise  Polynomial  Curve  Fitting  Algorithms 

In  this  section  we  present  our  algorithms  for  adaptively  finding 

smooth  piecewise  polynomial  uniform  approximations.  First,  we  define 

what  we  mean  by  a piecewise  polynomial  uniform  approximation.  Suppose 

the  function  to  be  approximated,  f,  is  defined  on  X a compact  set  of  real 

numbers,  and  [a,  b]  is  the  smallest  interval  containing  X.  A function  p 

is  a piecewise  polynomial  approximation  to  f of  degree  n if  there  are  m 

closed  intervals,  [a,  x.],  [x  , x„],  ...,  [xm  .,  b]=I  ,1  ,...,  I , such 

that  the  restriction  of  p to  I.,  is  in  = {p  : p is  a polynomial  of 

degree  less  than  or  equal  to  n}  for  each  j.  If  p restricted  to  1^  is  the 

best  uniform  approximation  to  f on  X 0 b from  IIn  for  each  j,  then  p is 

said  to  be  the  best  piecewise  polynomial  uniform  approximation  to  f of 

degree  n (with  respect  to  the  partition  1^,  I2,  ...,  1^  of  [a,  b]  0 X). 

See  [1]  and  [7]  for  a theoretical  study  of  this  problem  where  I,,..., I 

x m 

are  chosen  so  that  E,  = max  (|f(x)  - p(x)|}  is  the  same  for  all  j 

•'  xti.nx 

(thereby  giving  the  partition  containing  the  minimal  number  of  intervals 
needed  to  achieve  this  error). 

However,  in  many  applications  of  approximation  theory,  it  is  also 
desired  to  find  a smooth  approximation.  If  one  adds  this  constraint  to 
the  above  problem,  then  a best  smooth  piecewise  polynomial  uniform 
approximation  may  not  exist  having  the  same  error  on  each  subinterval. 
Even  if  such  an  approximation  does  exist,  the  number  of  subintervals  in 
the  partition  can  be  unreasonably  large.  Thus,  rather  than  finding  a 
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best  approximation  in  the  above  sense,  it  is  reasonable  to  first  stipu- 
late an  error  bound  for  ||  f - p||  and  a smoothness  constraint  on  the 
piecewise  polynomial.  Given  these  constraints,  one  then  should  seek  a 
partition  1^,  ...»  I of  [a,  b]  such  that  m is  in  some  sense  minimal  or 
small  and  the  piecewise  polynomial  uniform  approximation,  p,  to  f with 
respect  to  this  partition  satisfies  both  the  error  tolerance  and  smooth- 
ness demands. 

This  latter  approach  is  what  we  will  follow.  Typically,  our 
smoothness  constraint  will  be  that  p €.  CV[a,  b]  where  v = 0,  1 or  2. 

Our  adaptive  curve  fitting  algorithms  will  find  smooth  piecewise  poly- 
nomial uniform  approximations  satisfying  a prescribed  error  tolerance  and 
smoothness  requirements  by  adaptively  choosing  subintervals.  This 
partition  will  be  chosen  in  such  a way  that  the  error  of  the  best  uniform 
(or  restricted  range  approximation)  to  f on  the  intersection  of  each  sub- 
interval  with  X from  subject  to  interpolator^  constraints  guaranteeing 
the  desired  smoothness  is  less  than  the  prescribed  error  tolerance.  At 
the  same  time,  effort  is  made  to  keep  the  number  of  subintervals  as  small 
as  possible.  In  what  follows,  the  points  (x^}?_Q  will  be  called  knots, 
where  xQ  = a,  x^  = b,  and  1^  = [x^,  x^]. 


Algorithm  1 

We  assume  that  f fc  C(X).  We  use  the  following  notation:  n denotes 

the  degree  of  the  approximating  polynomials.  SMTH  denotes  the  number  of 
continuous  derivatives  desired  at  the  knots  (SMTH  < n);  TOL  denotes  the 
approximation  tolerance  we  wish  our  piecewise  polynomials  to  satisfy;  and 
LNGTH  denotes  an  approximate  minimum  length  we  will  allow  for  any 
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subinterval.  We  assume  these  values  are  all  given  and  fixed.  We  begin 

'Xj 

by  choosing  £ X to  be  the  largest  point  in  X which  satisfies 

(1)  x]_  - a >_  LNGTH 

(2)  If  p^  is  the  best  uniform  approximation  to  f from  on 

[a,  x ] f|  X,  then  sup  |f(x)  - p.(x)|  < TOL. 
xcta.^jnX 

If  x^  = b,  then  since  p^  is  a piecewise  polynomial  meeting  our  require- 
ments, the  algorithm  is  successfully  terminated.  If  no  such  x1  exists, 
then  the  algorithm  fails  and  is  terminated  with  an  appropriate  error 
message.  Otherwise,  we  choose  x^  to  be  the  largest  extreme  point  of 
f(x)  - p (x)  in  (a,  0 X (which  is  easily  available  to  us  from  the 

Remes  algorithm). 

We  continue  by  finding  successive  intervals  [x^  x2],  [x2,  x3],  ..., 

[x  b]>  and  polynomial  approximations  p0,  p.,  ....  p e n to  f so 

that  P^)(xi)  = P^^(xi)  for  i = 1,  2,  ....  m - 1,  j =0,  1,  ....  SMTH, 

and  so  that  sup  |f(x)  - p.(x)|  <_  TOL.  This  is  accomplished  as 

xe[xi_i,xi  ]0X 

follows . 

Suppose  we  have  determined  the  subintervals  [a,  x ],  [x^  Xj],  .... 
txi_2*  and  the  approximations  p^,  p2,  ...,  . Assume  further 

that  b - xj_^  _>  LNGTH.  We  now  determine  an  x.  and  a p^  meeting  the 
above  requirements.  We  begin  by  choosing  x.  c X to  be  the  largest  point 
in  X which  satisfies 

(1)  x£  - xi-X  > LNGTH, 

(2)  If  p.  is  the  best  approximation  to  f from  IIn  on  [x^  xj]  f)  x 

subject  to  the  constraint  that  p^(xi  ^ ), 


p 


1U 

j = o,  1,  ....  SMTH,  then  sup  ^ |f(x)  - p.(x.-|  <_  TOL. 

x£[x._i,xi](V<  1 

% 

If  x.  = b,  we  set  x^  = x^  = b,  and  the  algorithm  is  successfully  termi- 

nated.  If  no  such  x^  exists,  the  algorithm  fails  and  is  terminated  with 

an  appropriate  error  message.  Otherwise,  we  choose  x^  to  be  the  largest 

extreme  point  of  f(x)  - Pi(x)  in  (x.^  , 2L)  fl  X. 

Finally,  we  consider  the  special  case  where  b - x^  1 < LNGTH.  In 

this  case  we  change  x.  . to  some  x.  _ € X,  where  x.  _ < x.  , < x.  , < b. 

l-l  l-l  i-2  l-l  l-l 

So  that  b - x^_^  LNGTH.  Specifically,  we  choose  x^  to  be  a point  in 
X closest  to  (b  - x,.  ^ )/2  which  satisfies 

(1)  b - xi_1  >_  LNGTH, 

(2)  If  p^  is  the  best  approximation  to  f from  Hr  on  [x^  , b]  0 X 

subject  to  the  constraint  that  p^(x^  1>  = pp|(*,  ^ for 

j = 0,  1,  ....  SMTH,  then  sup  |f(x)  - Pi(x)|  _<  TOL. 

xrf£i,b]nx 

Again,  if  we  can  find  such  an  x^^  ^ , then  the  algorithm  is  successfully 
terminated.  If  not,  the  algorithm  fails  and  is  terminated  with  an 
appropriate  error  message. 


Algorithm  2 

We  assume  that  f,  l,  and  u € C(X)  and  satisfy  the  hypotheses 
necessary  for  restricted  range  approximation.  Let  n,  SMTH  and  LNGTH  be 
as  in  Algorithm  1.  If  TOL  is  the  tolerance  we  wish  the  approximation  to 
satisfy,  then  we  make  the  additional  assumption  that  1 and  u have  been 
chosen  so  that  if  p is  any  function  satisfying  i(x)  <_  p(x)  _<  u(x)  for  all 
x €.  X,  then  ||f  - p||  _<  TOL.  We  begin  by  choosing  x1  € X to  be  the 
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largest  point  in  X satisfying 

(1)  x1  - a > LNGTH, 

(2)  There  is  a best  restricted  range  approximation,  p^Cx),  to  f 

% 

from  nn  on  [a,  x^]  ft  X. 

'V. 

If  xx  = b,  then  the  algorithm  is  successfully  terminated.  If  no  such  x^^ 
exists,  the  algorithm  fails  and  is  terminated  with  an  appropriate  error 
message.  Otherwise,  we  choose  x^  to  be  the  largest  extreme  point  of 
f(x)  - p1(x)  in  (a,  x1>  H X. 

Proceeding  analogously  to  Algorithm  1,  in  general,  if  b - x. 
i LNGTH,  we  choose  x.  £ X to  be  the  largest  point  in  X so  that 

(1)  x.  - x.  > LNGTH, 

l i-l  — 

(2)  There  is  a best  restricted  range  approximation,  p.(x),  to  f 

from  nn  on  [x^ x^]  O X satisfying  the  constraints 
pij)(xi_l)  = j = 0,  1,  ....  SMTH. 

If  Xj  s b the  algorithm  is  successfully  terminated.  If  no  such  x. 
exists,  the  algorithm  fails  and  is  terminated  with  an  appropriate  error 
message.  Otherwise,  x^  is  chosen  to  be  the  largest  extreme  point  in 

(Xi-1’  *i^  x* 

If  b - x^_^  < LNGTH,  we  choose  x^ e X to  be  a point  closest  to 
(b  - x^  2 )/2  which  satisfies 

(1)  b - xi_  > LNGTH 

(2)  There  is  a best  restricted  range  approximation,  p.(x),  to  f 

(i  ) 

from  nn  on  [x^  pb]  ^ X satisfying  the  constraints  p^J'(>L 
= j = 0,  1 SMTH. 
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If  we  can  find  an  , the  algorithm  is  successfully  terminated;  if 
not,  the  algorithm  fails  and  an  appropriate  error  message  is  printed. 

Remarks 

We  first  note  that  in  algorithm  1 it  is  in  general  not  necessary 
to  compute  the  best  approximation  on  a particular  subinterval  in  order 
to  conclude  that  the  subinterval  is  too  long.  Indeed,  at  each  iteration 
of  the  Remes  algorithm  we  obtain  a lower  bound  for  the  error  of  the 
best  approximation  on  this  particular  subinterval,  namely  |a|. 
Consequently,  if  during  some  iteration  of  the  Remes  algorithm  |a|  > TOL, 
we  may  conclude  that  the  present  subinterval  is  too  large  and  terminate 
the  Remes  algorithm. 

In  our  implementation  of  these  algorithms,  the  x^  are  chosen  as 
follows.  At  each  step  of  this  iterative  procedure  we  will  let  ^ be  the 
current  largest  point  in  X such  that  requirements  (1)  and  (2)  of  the 
appropriate  algorithm  are  satisfied  on  [x^  a]  O X,  and  we  will  let  S' 
be  the  current  smallest  point  in  X such  that  fe  > a and  requirement  (2) 
of  the  appropriate  algorithm  fails  to  be  satisfied.  We  initialize  this 
process  by  computing  (or  attempting  to  compute)  the  appropriate  best 
approximation  on  [x^  , b]  D X.  If  this  approximation  satisfies 
requirement  (2),  then  we  set  x\  = b and  we  are  done.  If  the  approxi- 
mation  fails  to  satisfy  (2),  we  set  b = b.  Next,  let  t=min{x  € X : 
x - > LNGTH}.  If  the  approximation  on  [x.^,  t]  f\  X fails  to 

satisfy  (2)  then  the  algorithm  cannot  meet  the  desired  accuracy  and 
fails.  Otherwise,  we  set  a = t. 

In  general,  we  proceed  as  follows.  We  let  t = inf{x  €.  X : (^-a)/2 
< x < b).  If  this  set  is  empty,  we  set  t=sup{xC  X:a  <_  x <_  (b-a)/2r. 
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If  t = if  then  this  procedure  has  converged  and  we  set  x.  = t = a. 
Otherwise,  we  compute  (or  attempt  to  compute)  the  appropriate  approxi- 
mation on  [x.^,  t]  A X.  If  this  approximation  satisfies  (2)  then  we 
set  a = t.  If  this  approximation  fails  to  satisfy  (2)  then  we  set 
b = t.  We  continue  this  process  until  b - a is  less  than  some  prescribed 
tolerance,  at  which  point  we  accept  a as  a good  approximation  to  x.  and 
terminate  this  procedure.  We  compute  the  x^  in  a manner  analogous  to 
the  above. 

In  an  attempt  to  accelerate  the  convergence  of  the  above  scheme  in 
curve  fitting  algorithm  1,  we  have  tried  to  take  advantage  of  the  fact 
that  corresponding  to  each  a we  know  the  error  of  approximation  on 
[Xi  a]  fl  X,  call  it  SMLERR,  and  corresponding  to  each  b we  know  a 

'V. 

lower  bound  for  the  error  of  approximation  on  [x^  b] 0 X,  call  it 
BIGERR.  We  change  the  above  scheme  by  first  setting  a =(BIGERR  - TOL) 
/(BIGERR  - SMLERR)  and  then  setting  t = inf{xeX  : aa+(l-a)&<x<b}  or, 
if  this  set  is  empty,  we  set  t = sup(x  £X:a^_x<>o3'+(l-  S'Jb)  in 
the  general  iteration  described  above.  Hence,  if  SMLERR  is  very  close 
to  TOL,  t will  be  chosen  close  to  a.  Our  numerical  experience  has  shown 
that  this  procedure  only  works  well  when  approximating  uniformly  smooth 
functions  and  that  this  algorithm  cannot  be  significantly  improved  by 
allowing  the  Rentes  algorithm  to  run  to  completion  in  order  to  obtain 
the  true  error  of  approximation  for  BIGERR. 

If  f(x)  is  differentiable  on  some  interval  and  xQ  is  an  interior 
relative  extreme  point  of  f(x)  - p(x)  in  this  interval,  then  f’(xQ) 

= p'(Xq).  Hence,  by  "backing  off"  from  3f.  to  x^  we  hope  to  force  the 
first  derivative  of  the  approximating  polynomial  and  the  first 
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derivative  of  the  function  being  approximated  to  agree  closely  at  the 
knots,  thereby  dampening  oscillatory  effects  otherwise  prevalent  in  this 
type  of  approximation  when  we  require  the  approximation  to  be  at  least 
continuously  differentiable.  Our  numerical  experience  has  shown  that 
this  procedure  is  crucial  in  order  to  obtain  reasonable  approximations. 
If  we  only  want  our  approximations  to  be  continuous  at  the  knots,  then 
this  procedure  accomplishes  nothing  and  should  be  overridden. 

In  our  implementation  of  these  algorithms,  we  have  always  assumed 
that  X is  a finite  set  of  points.  Because  the  Remes  algorithm  can  only 
find  approximations  on  sets  of  at  least  n + 1-SMTH  points,  one  of  the 
uses  of  the  parameter  LNGTH  is  to  insure  that  every  set  on  which  we 
approximate  contains  at  least  n + 1-SMTH  points. 

If  we  set  fc(x)  = f(x)  - TOL  and  u(x)  = f(x)  + fOL  in  algorithm  2, 
the  approximations  obtained  from  algorithm  2 are  the  same  as  those 
obtained  from  algorithm  1.  However,  for  a broad  class  of  approximation 
problems,  the  additional  computational  complexity  of  algorithm  2 
justifies  the  use  of  the  simpler  algorithm  1.  We  will  discuss  this 
aspect  in  greater  detail  in  the  next  section.  Algorithm  2 is  particu- 
larly suited  to  problems  where  the  desired  tolerance  varies  over  the 
length  of  the  interval.  In  the  next  section  we  will  describe  a strategy 
which  essentially  allows  us  to  weight  data  points  which  contain  signi- 
ficant levels  of  noise  by  an  appropriate  choice  of  t(x)  and  u(x). 
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CHAPTER  III 

Numerical  Results 

Algorithms  1 and  2 have  been  implemented  as  FORTRAN  programs  and 
have  been  tested  on  Colorado  State  University's  CDC  CYBER  172.  Care 
has  been  taken  to  make  the  programs  reasonably  efficient  and  consistent 
with  each  other  so  that  CPU  time  comparisons  should  be  meaningful.  As 
examples,  the  functions  e^  on  [-1,  1],  |sin(x)|  on  [-n,  ir]  and  Jii 
on  [0,  2]  were  approximated  on  200  equally  spaced  points  with 
TOL  = .01,  SMTH  = 2,  and  N = 6.  When  using  algorithm  2 we  chose 
u(x)  = f(x)  + TOL  and  l(x)  - f(x)  - TOL.  Each  of  these  examples  is 
most  difficult  to  approximate  by  polynomials  near  x = 0.  Consequently, 
the  algorithms'  ability  to  automatically  decrease  the  length  of  the 
sub intervals  near  x = 0 and  then  to  recover  by  lengthening  them  for 
x > 0 is  tested.  As  noted  above,  the  approximations  computed  by  the 
two  algorithms  should  (and  do)  agree  up  to  machine  accuracy.  Below  is 
a table  listing  knot  locations  (subinterval  endpoints)  and  the  CPU  time 
in  seconds  used  by  each  algorithm  to  compute  the  piecewise  polynomial 
approximation. 
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,1*1 

|sin(x)  | 

<&eS> 

Algorithm  1 

.979 

1.095 

.876 

Algorithm  2 

1.518 

1.539 

1.367 

-1.0 

-3.14 

6.0 

- .347 

- .268 

.0804 

- .0653 

.0474 

.231 

Knot 

.0452 

.205 

.352 

Locat ions 

.196 

.458 

.814 

.437 

.868 

2.0 

1.0 

1.53 

3.14 

We  have  found  that  in  order  for  these  algorithms  to  be  stable,  it 
is  crucial  that  when  computing  smooth  approximations  (SMTH  > 0),  the 
derivative  of  the  approximating  polynomial  agree  closely  with  the  slope 
(or  approximate  derivative)  of  the  function  being  approximated  at  the 
right  endpoint  of  every  subinterval  except  the  last.  When  approximating 
on  a dense  enough  subset  of  an  interval,  the  strategy  of  choosing  x^  to 
be  the  last  interior  extreme  point  of  f - p given  in  the  description  of 
the  algorithms  seems  to  be  sufficient.  However,  for  more  widely  spaced 
points,  the  true  last  interior  extreme  point  of  f - p,  x*,  generally  is 
not  in  X and  often,  the  closest  point  in  X to  x*,  which  is  found  by  the 
Remes  algorithm,  is  not  a point  at  which  p'  agrees  very  closely  with  the 
slope  of  f.  Our  numerical  experience  has  shown  that  by  setting  f'(x) 
to  be  the  derivative  of  the  centered  quadratic  interpolation  of  f at  x 
and  choosing  x^  to  be  the  largest  extreme  point  returned  by  Remes  such 
that  (f'(x)  - p'(x)}  is  a minimum  we  obtain  a considerably  more  stable 
algorithm. 

For  some  applications  it  might  be  preferable  to  be  able  to  choose 
ahead  of  time  a particular  knot  location  and  to  specify  interpolatory 
constraints  at  these  knots.  For  example,  when  approximating  |sin(x)  | 
it  might  be  advantageous  to  force  a knot  to  be  located  at  x = 0 and  to 
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require  p(0)  = 0,  p'(0)  = -1  as  x approaches  0 from  the  left,  and 
p'(0)  = -1  as  x approaches  0 from  the  right.  The  modifications  to  the 
algorithms  as  given  necessary  to  accomplish  this  would  not  be  extensive 
or  difficult. 

Because  uniform  approximations  weight  each  data  point  equally,  they 
are  particularly  suited  for  approximating  precise  mathematical  functions 
or  for  approximating  data  with  noise  levels  that  are  very  small  relative 
to  the  desired  accuracy.  However,  when  these  algorithms  are  used  to 
approximate  functions  containing  considerable  levels  of  noise,  the 
approximations  tend  to  follow  the  noise  patterns  more  than  is  desirable 
so  that  near  abrupt  changes  in  the  data,  the  approximations  tend  to 
begin  to  oscillate  and  generally  it  requires  several  subintervals  to 
dampen  these  oscillations. 

For  example,  we  were  given  some  experiemtnal  data  involving  the 
release  of  bitumen  from  oil  shale  heated  to  a constant  temperature  as  a 
function  of  time.  Because  relatively  lew  data  points  were  available, 
we  filled  in  the  gaps  between  the  data  points  by  linear  interpolation  so 
that  we  approximated  on  200  equally  spaced  points.  Figure  1 is  a plot 
of  the  data  being  approximated,  the  constraining  curves  u(x)  and  *(x), 
and  the  approximation  generated  by  algorithm  2 where  N = 6,  SMTH  =2, 

TOL  = 2.5,  u(x)  = f (x)  + TOL,  and  i(x)  = f(x)  - TOL.  Notice  that  as 
oscillations  appear  in  the  data  between  times  of  22  and  34  minutes, 
greatly  exaggerated  oscillations  appear  in  the  approximation. 
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RESTRICTED  RANGE  75  GAL/TON  TEMP=425  BITUMEN 
N=6,  NSMTH=2,  T0L=2.500.  Knots  are  indicated  by  0. 


Figure  1 

In  order  to  dampen  this  effect,  we  have  found  that  by  appropriately 
choosing  u(x)  and  j,(x)  we  can  improve  upon  this  approximation. 
Essentially,  we  choose  u(x)  and  £(x)  such  that  for  each  x 6.  X, 
max{u(x)  - f ( x ) , f(x)  - *.(x)}  = TOL,  and  u(x)  and  l(x)  are  as  smooth 
as  is  reasonably  possible.  Thus,  the  data  points  are  scattered  inside 


a band  of  varying  width,  but  neither  boundary  is  further  than  TOL  from 
a data  point.  When  one  data  point  differs  sharply  from  those  on  either 
side  of  it,  we  appropriately  choose  either  u(x)  or  t(x)  to  be  close  to 
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this  data  point.  The  effect  of  this  choice  of  u and  l is  to  tend  to 
force  the  approximation  toward  the  more  "believable"  side  of  the  data 
which  results  in  fewer  oscillatory  problems.  Figure  2 is  a plot  of  the 
same  example  shown  in  Figure  1,  but  with  u(x)  and  H( x ) chosen  as  above. 

RESTRICTED  RANGE  75  GAL/TON  TEMP=425  BITUMEN 
N=6,  IJSMTH=2 , T0L=2 . 500 . Knots  are  indicated  by  0. 


Figure  2 

If  a technician  collecting  data  knows  the  general  trend  his  data 
should  follow,  he  can  set  u(x)  and  £(x)  accordingly,  forcing  the 
approximation  to  lie  within  theoretical  bounds  instead  of  just 
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prescribing  an  error  tolerance.  Or,  if  he  knows  that  part  of  his  data 
is  highly  accurate  relative  to  the  rest,  he  can  easily  vary  the  toler- 
ance he  is  willing  to  accept  accordingly.  If  he  finds  that  the  approxi- 
mation generated  using  one  set  of  restraining  curves  allows  the  approxi- 
mation to  have  undesirable  characteristics,  generally  he  can  alter 
these  characteristics  by  appropriately  changing  u(x)  and  £(x). 

There  is  no  requirement  inherent  in  these  algorithms  for  the  data 
points  to  be  equally  spaced.  By  adding  more  points  in  regions  where 
there  are  quick  variations  in  the  function  being  approximated,  smoother 
approximations  result  due  to  the  fact  that  the  algorithm  will  find  more 

optimal  right  end  points  as  they  "back  off"  from  the  x.. 

2 . . 

Because  L approximations  minimize  the  effects  of  normally  distri- 
buted random  noise,  we  are  experimenting  with  altering  algorithm  1 by 
. 2 

computing  best  L approximations  in  each  subinterval  instead  of  best 

2 

uniform  approximations,  the  theory  for  best  L restricted  range 

approximations  is  very  difficult  and  has  not  been  developed;  therefore, 

2 

we  cannot  construct  the  analog  to  algorithm  2 using  best  L approxima- 
tions. We  will  report  on  this  work  in  a future  paper. 

Recently,  John  Rice  has  written  several  papers  in  which  he  describes 
an  adaptive  piecewise  polynomial  algorithm  which  uses  local  Hermite 
interpolation  instead  of  uniform  approximation  on  each  subinterval  to 
obtain  each  polynomial  piece  (see  for  example  [9]).  His  algorithm 
requires  that  the  first  SMTH  derivatives  or  approximate  derivatives  of  f 
be  available  in  order  to  compute  approximations  which  have  SMTH  contin- 
uous derivatives.  Rice's  adaptive  strategy  also  differs  from  ours;  he 
uses  a bisection  strategy  which  can  be  described  as  follows.  First, 


compute  the  Hermite  interpolating  polynomial  to  f on  [a,  b].  That  is, 
compute  the  polynomial  which  interpolates  f and  its  first  SMTH  deri- 
vatives at  both  endpoints,  as  well  as  interpolating  f at  k evenly 
distributed  points  in  (a,  b),  where  k = DEGREE  - 2* SMTH  - 1,  and  DEGREE 
is  the  degree  of  the  approximating  polynomial.  Next,  measure  the  error 
of  approximation  using  any  preselected  iJ3  norm  (p  >_  1).  If  the  error 
is  less  than  TOL,  the  algorithm  terminates,  otherwise  bisect  [a,  b] 
into  two  subintervals  and  check  to  see  if  the  Hermite  interpolating 
polynomial  on  [a,  (a  + b)/2]  differs  from  f in  norm  by  no  more  than  TOL. 
Meanwhile,  place  the  subinterval  [(at  b)/2,  b]  in  a "stack"  to  be 
processed  later.  Continue  bisecting  subintervals,  placing  the  right 
half  of  the  interval  on  the  top  of  the  stack  to  be  processed  later  and 
check  the  left  half  to  see  if  the  error  of  approximation  by  the  Hermite 
interpolating  polynomial  is  less  than  TOL.  When  a subinterval  and  its 
associated  approximation  are  found  which  meet  the  desired  tolerance, 
we  accept  this  approximation  as  a "piece"  of  the  piecewise  polynomial 
approximation  and  continue  the  algorithm  by  removing  the  subinterval 
which  is  currently  at  the  top  of  the  stack  and  repeating  the  above  on 
it  unless  the  stack  is  empty,  at  which  point  terminate  the  algorithm. 

Rice's  routine  requires  the  value  of  the  function  and  its  deriva- 
tives at  very  many  points  in  [a,  b],  even  though  it  may  not  actually 
access  these  values.  His  algorithm  is  particularly  suited,  then,  for 
approximating  precise  mathematical  functions  for  which  this  data  is 
readily  available.  Our  algorithms  seem  more  suitable  for  data  reduction 
in  the  sense  that  instead  of  storing  a large  number  of  data  points, 
the  user  could  instead  store  the  coefficients  of  several  polynomials 
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resulting  in  a savings  of  storage  space.  Our  approximations  seem 
preferable  for  this  purpose  since  our  algorithms  do  not  require  any 
information  about  derivatives,  and  no  special  configuration  of  the  data 
is  needed  for  our  adaptive  strategy  to  work. 

By  using  interpolatory  constraints  on  both  endpoints  of  each 
subinterval,  and  interpolating  the  values  of  the  derivative  of  f at 
these  points,  we  could  modify  our  algorithms  to  use  a bisection  adaptive 
strategy.  However,  in  order  to  take  greatest  advantage  of  our  use  of 
the  best  uniform  approximation  operator  as  opposed  to  the  Hermite 
interpolation  operator,  we  felt  that  it  was  preferable  to  use  the 
majority  of  our  degrees  of  freedom  for  inside  subintervals  instead  of 
using  them  up  on  interpolation  requirements  at  the  endpoints.  Essen- 
tially, in  our  algorithms  we  traded  a somewhat  faster  run  time  (since 
the  uniform  approximation  operator  is  computationally  very  slow  relative 
to  the  Hermite  interpolation  operator)  for  fewer  knots,  the  freedom  from 
needing  to  know  derivative  information,  and  in  algorithm  2,  greater 
control  of  the  characteristics  of  the  approximations. 
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